Photon statistics without counting photons 
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We show how to obtain the photon distribution of a single-mode field using only avalanche pho- 
todetectors. The method is based on measuring the field at different quantum efficiencies and 
then inferring the photon distribution by maximum-likelihood estimation. The convergence of the 
method and its robustness against fluctuations are illustrated by means of numerically simulated 
experiments. 
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INTRODUCTION 



Optical signals and nonclassical states of light have 
been the subject of constant attention over the last three 
decades. The photon distribution, besides fundamental 
interest, plays a major role in quantum communication 
schemes based on light beams, and this stimulated mam/ 
experiments on the statistical properties of radiation [lj . 
More recently, the rapid development of quantum infor- 
mation processing once again posed the issue of effec- 
tive methods to investigate photon statistics 0, 01 ■ Pho- 
ton distribution is usually obtained by photon counting 
in time intervals much shorter than the coherence time 
of the beam under investigation. A long time stabil- 
ity is thus needed and this requirement becomes more 
and more strict as far as the intensity of light becomes 
lower. As a matter of fact, most photon counting exper- 
iments involved laser beams. Effective photon counters 
have been also developed, though their current operating 
conditions are still extreme 0. 

The advent of quantum tomography provided a novel 
method to measure photon distribution y. However, 
the tomography of a state, which have been applied to 
several quantum states [7} , needs the implementation of 
homodyne detection, which in turn requires the appro- 
priate mode matching of the signal with a suitable local 
oscillator at a beam splitter. 

In this paper, we address a simple method to obtain 
the photon distribution without directly counting pho- 
tons. In this scheme, repeated preparations of the signal 
are revealed through avalanche photodetectors (APD) 
at different quantum efficiencies. The resulting on/off 
statistics is then used to reconstruct the photon distri- 
bution through maximum-likelihood estimation. Since 
the model is linear and the photon distribution is a set 
of positive numbers, then the maximum of the likelihood 
functional can be found iteratively by the expectation- 
maximization (EM) algorithm || Q . The method does 
not require long time stability and involves only simple 
optical components. The number of experimental runs 



depends on the signal under investigation, roughly in- 
creasing with its nonclassicality. 

The idea of inferring photon distribution through de- 
tection at different efficiencies has been already ana- 
lyzed theoretically |lfj |. and implemented to realize a 
multichannel fiber loop detector ^l|- Here we describe 
other possible implementations, analyze the reconstruc- 
tion when only a subset of values < ?y m i n < t] < r) max < 
1 of the quantum efficiency is available, and discuss in 
details the statistical properties of the method: conver- 
gence and robustness against fluctuations. 

The paper is structured as follows: In Section[H]we in- 
troduce the problem, and show that simple estimation by 
inversion cannot be used due to lack of precision. Possi- 
ble implementations of the measurement scheme are also 
described. Then, in Section UTTI we illustrate the recon- 
struction of photon distribution by iterative solution of 
maximum likelihood estimation. In Section llVI the con- 
vergence properties of the method, as well as its robust- 
ness to fluctuations, are discussed on the basis of several 
Monte Carlo simulated experiments, performed on differ- 
ent kind of signals. Section[V]closes the paper with some 
concluding remarks. 



II. ESTIMATION BY INVERSION 

Given a single-mode state g — J2 n m Qnm\n)(m\ we are 
interested in the photon distribution, i.e in the set of 
positive numbers g n = g nn > 0. We assume to have at 
disposal APDs, which can only discriminate the vacuum 
from the presence of radiation, with a certain quantum 
efficiency. This kind of measurement, on/off photodetec- 
tion, is described by the following probability operator- 
valued measure (POVM) {n o ff,n on } 

oc 

n o ff(ry) = 5>-r?)>Xn| 

n=0 

n oa (ry) = E-n off (1) 
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77 being the detector's quantum efficiency, i.e. the prob- 
ability that an incoming photon lead to a click of the 
detector. For any given state the detector does not click 
with a probability p s{v) = Tr{pn o ff (77)}, that reads as 



2 



follows 



Pos(v) = ^0--v) n Qv 

n=0 



(2) 



From now on we suppress the subscript "off" and always 
mean p a s when we write p. The "off" probabilities for 
a set of N detectors measuring the same quantum state 
with different quantum efficiencies are then 



Po0?o) = E„(l - Vo)" Qn 

PlM = J2n( 1 -Vl) n Qn 
Pn(Vn) = £„(! - Vn)" Qn 



(3) 



If we know all of the r\ v 's values, equation (|3J) is a linear 
system with unknowns {g n }- In practice, it is not nec- 
essary to have at disposal many detectors with different 
quantum efficiencies, since a suitable tuning of rj can be 
obtained by optical filters or through an interferometric 
setup. In fact, besides the fiber loop scheme of ^jj, dif- 
ferent quantum efficiencies can be obtained inserting a 
set of optical filters before the detector, or by the scheme 
of Fig. ^ where a single APD is needed, and lower effi- 
ciencies are obtained by varying the internal phase-shift <p 
of the interferometer. Since the overall transmissivity of 
the interferometer is r = cos 2 </>, and a tuning of <j> of the 
order of 7r/500 can be actually achieved, we may safely 
assume that about 100 different values of r\ = between 
?y m i n ~ and T] max = t/apd can be obtained. 



and the coefficients matrix V (for r\i ^ r\j Vi,j) is a 
nonsingular Vandermonde matrix of order n. If we put 
Xi = 1 — rji the V matrix reads 

,2 _rT— 1 



V = 



X 
Xi 



2 

-1 XjT- 



(7) 



and the photon distribution can be obtained by matrix 
inversion g = V -1 p. The same approach can be used 
also when we are able to reveal the state with N > n 
number of ry's. In this case the system in Eq. should 
be solved in the least square sense leading to g = W p, 
where W = (V T V) _1 V T is the Moore-Penrose inverse of 
V. 

Unfortunately, the reconstruction of g n by matrix in- 
version cannot be used in practice since it would require 
an unreasonable number of experimental runs. In fact, 
most of the quantities Xi entering the expression of V are 
of order 10 _1 , and therefore the V _1 's entries in the j-th 
line are of order a; - ^ -1 ' . This means that the reconstruc- 
tion of Qj requires the multiplication of the experimental 
frequencies Pi by quantities of the order x~^~^', which 
in turn implies that a sound reconstruction of gj needs 
that pi must be precise at least to the (j — l)-th deci- 
mal digit, i.e. a minimum of lO^ -1 experimental runs. 
In addition, for increasing ./V and n the inversion of 
must be done numerically leading to errors that quickly 
become unacceptably large. 




FIG. 1: A possible experimental setup for simulating different 
quantum efficiencies. The signal passes through an interfer- 
ometer with internal phase-shift (j) and is then revealed by a 
high-efficiency APD. The transmissivity of the interferometer 
is r = cos 2 (j>, and the overall efficiency of photodetection is 

V = T77APD- 



Suppose now that the g n 's are negligible for n > n 
and that we are able to measure the signal with N = n 
different 77's. In this case equation © is a linear system 
of the form 



• Q 



where 



P 
Q 



{P0,Pl, ■ ■ ■ ,Pn-l} 
{g , Ql,---, gn-l} 1 



(4) 

(5) 
(6) 



III. MAXIMUM-LIKELIHOOD ESTIMATION 

The problems illustrated in the previous Section can be 
circumvented by considering equation @ as a statistical 
model for the parameters g n to be solved by maximum- 
likelihood (ML) estimation. We assume N > n and, for 
sake of simplicity, we define 



Pu = Pu{r}u), 
so that equations J3Jl can be rewritten as 



Pu 



. A-wn Qn 



(8) 
(9) 



Since the model is linear and the parameters to be esti- 
mates are positive (LINPOS problem), then the solution 
can be obtained using the Expectation-Maximization 
algorithm (EM) 0, By imposing the restriction 

En Bri = 1, we obtain the iterative solution 



qI: +1) - $ E 



A h 



(10) 



where Qn is the value of g n evaluated at i-th iteration, 
C„ = n v ^2 m A um , n„ being the total number of exper- 
imental runs with rj — rj v , \\ v in the number of no- click 
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events for rj — rj v and p u [{Qn }] are the frequencies p v 
calculated using the reconstructed distribution {g„ } at 
the «-th iteration. By introducing the symbol j v = \\ v jn v 

for the experimental frequencies, the expression of 
rewrites as 



L^™ P„[{0n }] 



(11) 



EM algorithm is known to converge unbiasedly to the 
ML solution. Indeed, it has been already used to infer 
the photon distribution from random phase homodyne 
data |12|. The confidence interval on the determination 
of the element g n can be given in terms of the variance 
cr„ = 1/ y/N F n , N being the number of measurements 
and F„ the Fisher's information 113 



where 



<lv 



(12) 



are the renormalized probabilities of no-click with quan- 
tum efficiency r\ v . Nq — Ylun A-vnQn is the global fraction 
of no-click events (irrespective of the quantum efficiency). 

Notice that Eq. (| 1 1 11 provides a solution once an ini- 
tial distribution {f?n°^} is chosen. In our simulated ex- 
periments we start from the uniform distribution = 
(1 + n) _1 in [0,n]. Other choices, the only constraint 

being Qn ^ =/= 0, Wi, do not dramatically influence the 
convergence properties of the algorithm. 



IV. MONTE CARLO SIMULATED 
EXPERIMENTS AND DISCUSSION 



from the actual probabilities as calculated from the the- 
oretical photon distribution. As a measure of accuracy 
we adopt the fidelity 



Qn Q 



(k) 



(15) 



between the reconstructed distribution and the theoret- 
ical one. In Figs. 121 II (right) we report versus the 
number of iterations for different signals. As it is ap- 
parent from the plots, the total error is a good marker 
for the convergence of the algorithm, while the normal- 
ization factor — J2 n £™ ^ 1 (ideally zero at each 
step) is not. Notice, however, that the minimum of the 
total error does not always coincide with the maximum 
fidelity of reconstruction (Figs. El and 01 , which means 
that our method is slightly biased, especially for fast re- 
construction, i. e. when it converges quickly as it happens 
for coherent signals. We have numerically observed that 
this problem can be circumvented by using a number of 
iterations n;t — n x approximately equal to the number 
of data. Currently we have no precise explanation of this 
phenomenon, and provide it as a heuristic prescription 
leading to best performances for a large class of quan- 
tum signals. 
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We have performed several numerical simulations in 
order to check the accuracy and reliability of our method 
by varying the different parameters. Since our solution of 
the ML estimation is obtained from an iterative solution, 
the most important aspect to keep under control is its 
convergence. As a measure of convergence we use the 
total absolute error at the fc-th iteration 

N 

£ «=5>W|, (13) 

v=0 

where 

n-l 

4 fe) = Pv PvliQ^}} =Pv- £(i - wTq^ ■ (14) 

The total error measures the distance of the prob- 
abilities Pv[{Qri^}], as calculated at the fc-th iteration, 



FIG. 2: Reconstruction of the photon distribution of a coher- 
ent state | a) with average number of photons (eJa) = \a\ 2 — 
5.20. On the left: the reconstructed photon distribution af- 
ter nit = 10 4 iterations. On the right: the normalization 
factor S^ k ' and total error e^ k ' il l as a function of number 
of iterations. The confidence interval has been evaluated as 
l/V n xF n where F n is the Fisher information. The number 
of simulated data and the number of iterations are given by 
fix = Wit = 10 5 . In all the simulated experiments N = 50 dif- 
ferent quantum efficiencies have been used with a minimum 
efficiency ry m i n = 0.02. The maximum efficiency is given by: 
(a) 7| m „ = 0.99; (b) 7/ m ax = 0.5. The Hilbert space is trun- 
cated at n = 20. 

We have performed simulated experiments for coher- 
ent states | a) = D(a)\0), squeezed states |a, £) = 
D(a)S(£)\Q) and superposition of Fock states \ip) = 
2~ 1 / 2 (|n 1 ) + |n 2 )), where D{a) = exp(ao 1 ' — aa) is the 
displacement operator and S(£) = exp(i£a^ 2 — \£,a 2 ) 
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is the squeezing operator. Squeezed states have been 
parametrized through the total average photon number 
and the squeezing fraction 



a* a) 

c 



+ I£l7(i 

|a|7<a+a) 



lei 2 ) 



(16) 



The value £ = 1 corresponds to a squeezed vacuum and 
C = to a coherent state. As shown in Fig. [3 the al- 
gorithm converges quite fast for coherent states, while 
for nonclassical states such as squeezed states (Fig. |3J) 
the number of needed iterations is larger. The right 
plots in Fig. also indicate that increasing N does 
not always improve accuracy. In Fig. 0] we show re- 
construction for the unbalanced superpositions of Fock 
states \i> 2 ) = (2/3) 1 /2| 2 ) + {l/Z) 1 / 2 ^). 

Concerning the values of the quantum efficiency, we 
used N values of 77 uniformly distributed in [?7 m i n , r) max ] 
with r] m i n ~ and f? max < 1. In principle, a different dis- 
tribution (not uniform) may influence the performances 
of the algorithm. We found, however, that both conver- 
gence and accuracy are not much affected by a different 
choice, which may become relevant only if the spacing 
between the efficiency values becomes smaller. It should 
be noticed that the algorithm works well also when ?7 max 
is considerably smaller than unit. This is a relevant fea- 
ture of the method in view of its experimental imple- 
mentations in different working regimes. In Figs. 121 II 
(bottom left) we report the reconstructions obtained as- 
suming ?7 max = 0.5 (coherent states and superpositions) 
and r/ max = 0.7 (squeezed states). 
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FIG. 4: Reconstruction of the photon distribution of a unbal- 
anced superpositions of number states (see text) . The number 
of simulated data and the number of iterations are given by: 
(a) n x = 10 4 ,nit = 10 6 ; The maximum efficiency is given by: 
(a) 7/ ma x = 0.99; (b) 7? max = 0.5. The other experimental 
parameters are the same as in Fig. |2] 



case we have checked that the algorithm is able to recon- 
struct accurately the norm of the included part, such that 
a simple check of the distribution norm allows to optimize 
n (and in turn N) in few steps. This is a remarkable fea- 
ture of the algorithm, since in general a large n improves 
convergence but doesn't guarantee better accuracy. 




FIG. 5: Fidelity G'*' versus the number of iterations. Left: 
for a squeezed state with (a T a) = 1.0 and different squeezing 
fractions £ (left). Right: for a squeezed state with (a? a) = 1.0, 
£ = 0.75 and different numbers N of 77's values. In both cases 
the maximum number of iterations is n-,t = 10 6 . 



FIG. 3: Reconstruction of the photon distribution of a 
squeezed state with squeezed fraction £ = 0.99 and average 
photon number (a'a) = 0.5. Notice that a larger number of 
iterations is needed in comparison with the coherent signals' 
case. The number of simulated data and the number of iter- 
ations are given by n x — 10 5 ,riit = 5 x 10 . The maximum 
efficiency is given by: (a) ?) max = 0.99; (b) ?7 max = 0.7. The 
other parameters are the same as in Fig. [5] 

In experiments where we have no a priori information 
on the state under investigation it could happen that 
part, or even most, of the number distribution g n lies 
outside the reconstruction region (from to n). In this 



The error bars in the plots have been calculated using 
the Fisher information (|12f) . that explicitly reads: 

where A^o is the total number of no click events. 

A question may arise about the robustness of the 
method against fluctuations in the value of the r\ v (which, 
in the case of the interferometric implementation of Fig. 
^ may occur as a consequence of phase fluctuations), 
i.e. whether or not their precise knowledge is needed. 
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FIG. 6: Fidelity versus the number of iterations for a 

squeezed state with {a^a) = 1.5 and £ = 0.75. Each line 
represent a different simulated run. Left: n x = 10 5 , right: 
n x = 10 6 . In both the plot the maximum iterations number 



is rut = 10° 
Fig.H 



while the other parameters are the same as in 



In order to check robustness we have performed simu- 
lated experiments where, during the run, the quantum 
efficiency may fluctuate. In particular, we assumed each 
T\ v uniformly distributed in the range (—a + a + rj v ), 
where 



Vn 



aN 



(18) 



and a is a positive number. The value a = 2 corresponds 
to each rj^ fluctuating in an interval as large as the spac- 
ing {r]u+i — i]v) around its expected value. The values 
of p v change accordingly during the run. Our results are 
summarized in Fig. [7| The reconstruction is not dra- 
matically affected by fluctuations, though errors bars are 
slightly larger. We conclude that the method is robust 
against fluctuations. 



CONCLUSIONS 



simple optical components, and allows reconstruction 
with APD quantum efficiency considerably smaller than 
unit. The convergence of the method, and its robust- 
ness against fluctuations of quantum efficiency have been 
demonstrated numerically, by means of Monte Carlo sim- 
ulated experiments. 
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FIG. 7: Reconstruction of the photon distribution for different 
signals for fluctuating r\ v . (a, a'): coherent state with (a* a) = 
5.20; (b,b'): squeezed state with (o t a> = 0.50 and £ = 0.99. 
The number of simulated data and the number of iterations 
are given by: (a,&')n x = 10 5 ,nj t = 10 5 ; (b,b') n x = 10 6 ,riit = 
5 x 10 6 ,(a* a)fl t = 5.01. In all the simulated experiments we 
set a — 2 in Eq. 1181 . The maximum efficiency is given by: 
(a,b) rj 

iii • :. = 0.99; (a') r\ max = 0.5; (b') rj max 

= 0.7. The other 
experimental parameters are the same as in Fig. [5] 



We analyzed in details an iterative algorithm to in- 
fer the photon distribution of a single-mode radiation 
field using only avalanche photodetectors. The method 
is accurate and statistically reliable for a large class 
of Gaussian (coherent and squeezed) and non Gaussian 
states (superpositions and mixtures of |n) states), pro- 
vided that on/off photodetection may be performed at 
different quantum efficiencies. The scheme involves only 
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